Anàlisi de dissenys factorials de dos o més factors creuats
Antonio Miñarro (aminarro@ub.edu)
12 novembre, 2024
1 Anàlisi de dissenys factorials de dos o més factors creuats
Podeu accedir a una versió electrònica actualitzada d’aquest document clicant al següent enllaç aminarroa.github.io/DEAD.html
1.1 Anova de dos factors fixos
1.1.1 Exemple 1: Analitzant dues causes de variabilitat
Un investigador vol avaluar l’eficàcia de 4 possibles vacunes per a una malaltia. L’investigador opina també que el número de dosis subministrades: una, dos o tres, pot tenir una influència en el nivell de protecció aconseguit. Efectuant un estudi en rates de laboratori opta per dividir un total de 48 rates en grups de 4 i sotmetre a cada grup a una vacuna i un número de dosis prefixades. Es mesura la quantitat d’anticòs Ig G en gr/l present a l’organisme.
anticos<-c(6.1,8.2,6.3,6.5,7.5,11.0,9.5,7.1,6.6,8.8,6.3,6.6,4.3,5.2,5.6,6.2,
3.6,5.2,4.4,5.6,2.9,5.1,3.5,10.2,4.0,4.9,3.1,7.1,2.3,12.4,4.0,3.8,
2.2,3.0,2.3,3.0,2.1,3.7,2.5,3.6,1.8,3.8,2.4,3.1,2.3,2.9,2.2,3.3)
vacuna<-gl(4,12)
dosis<-factor(rep(1:3,times=16))
immuno<-data.frame(anticos,vacuna,dosis)L’objectiu és estudiar simultàniament si els dos factors vacuna i dosi, amb 4 i 3 nivells respectivament, poden influir en la variable resposta quantitativa quantitat d’anticòs.
Els individus han estat aleatoritzats respecte els factors.
1.1.2 Disseny de dos factors amb interacció. Estructura de les dades
L’estructura de les dades correspon a una taula bidimensional on cada cel·la equival a una combinació de nivells dels dos factors, i dintre de cada cel·la tenim les rèpliques que corresponen a la condició experimental que representa la combinació.
\[ \begin{array}{cccc} & B_1 & \cdots & B_b \\ A_1 & y_{111},...,y_{11n_{11}} & \cdots & y_{1b1},...,y_{1bn_{1b}}\\ \vdots & \vdots & & \vdots \\ A_a & y_{a11},...,y_{a1n_{a1}} & \cdots & y_{ab1},...,y_{abn_{ab}} \end{array} \]
Si el disseny és balancejat \(n_{ij} = n\) per tot \(i=1,...,a\) i \(j=1,...,b\).
D’ara endavant si no diem el contrari assumirem que el disseny és balancejat.
1.1.3 Model lineal de l’ANOVA de dos factors fixos creuats
El model lineal que respon a l’enunciat de l’experiment és:
\[ Y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]
on
\(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),
\(\mu\) és la mitjana general,
\(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A,
\(\beta_j\) és un paràmetre del model que mesura l’efecte del nivell \(j\) del factor B,
\(\gamma_{ij}\) és un paràmetre del model que mesura l’efecte de la interacció del nivell \(i\) del factor A amb el nivell \(j\) del factor B (també és freqüent representar-lo com \((\alpha \beta)_{ij})\),
\(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(E(\epsilon_{ijk}) = 0\) i \(Var(\epsilon_{ijk})=\sigma^2\).
Assumirem també les següents restriccions sobre els paràmetres \(\sum_{i=1}^a \alpha_i = 0\), \(\sum_{j=1}^b \beta_j = 0\), \(\sum_{i=1}^a \gamma_{ij} = 0 \,\, j=1,...,b\) i finalment \(\sum_{j=1}^b \gamma_{ij} = 0 \,\, i=1,...,a\).
Com veurem serà necessari també assumir una distribució normal sobre els errors
\[ \epsilon_{ijk} \sim N(0,\sigma)\]
i per tant cada observació experimental
\[ Y_{ijk} \sim N(\mu_{ij},\sigma) \]
on \(\mu_{ij} = \mu + \alpha_i + \beta_j + \gamma_{ij}\).
La interacció la podem interpretar com l’efecte de la combinació de dos nivells, un de cada factor, que provoca una modificació en la variable resposta que va més enllà de l’efecte merament additiu.
1.1.4 Contrast d’hipòtesis
L’Anova de dos factors creuats amb interacció permet contrastar les següents hipòtesis:
- Hipòtesi sobre la significació del factor A
\[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]
- Hipòtesi sobre la significació del factor B
\[ H_0: \beta_1 = \beta_2 = \cdots = \beta_b = 0 \] \[ H_1: \beta_j \ne \beta_{j'} \mbox{ per algun } j \ne j' \,\,\, ( \mbox{ equivalentment } H_1: \exists \beta_j \ne 0)\]
- Hipòtesi sobre la significació de la interacció
\[ H_0: \gamma_{11} = \gamma_{12} = \cdots = \gamma_{ab} = 0 \] \[ H_1: \exists \gamma_{ij} \ne 0\]
1.1.5 Estimacions dels paràmetres del model
1.1.5.1 Definicions
– \(N=abn \,\,\,\) (Disseny balancejat)
– \(Y_{i..} = \sum_{j=1}^{b} \sum_{k=1}^{n}Y_{ijk}\)
– \(Y_{.j.} = \sum_{i=1}^{a} \sum_{k=1}^{n}Y_{ijk}\)
– \(Y_{ij.} = \sum_{k=1}^{n} Y_{ijk}\)
– \(Y_{...} = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} Y_{ijk}\)
– \(\overline{Y}_{i..}= \frac{Y_{i..}}{bn}\)
– \(\overline{Y}_{.j.}= \frac{Y_{.j.}}{an}\)
– \(\overline{Y}_{ij.}= \frac{Y_{ij.}}{n}\)
– \(\overline{Y}_{...}= \frac{Y_{...}}{N}\)
1.1.5.2 Estimadors puntuals dels paràmetres del model
– \(\hat{\mu} = \overline{Y}_{...}\)
– \(\hat{\alpha}_i = \overline{Y}_{i..} - \overline{Y}_{...} \,\,\,\, (\alpha_i = \mu_i - \mu)\)
– \(\hat{\beta}_j = \overline{Y}_{.j.} - \overline{Y}_{...} \,\,\,\, (\beta_j = \mu_j - \mu)\)
– \(\hat{\gamma}_{ij} = \overline{Y}_{ij.} - \overline{Y}_{i..}- \overline{Y}_{.j.}+ \overline{Y}_{...} \,\,\,\, (\gamma_{ij} = \mu_{ij} -\mu_i -\mu_j + \mu)\)
1.1.5.3 Valor predit i residu
– \(\hat{Y}_{ijk} = \overline{Y_{ij.}}\)
– \(\hat{e}_{ijk} = Y_{ijk} - \overline{Y_{ij.}}\)
1.1.6 Descomposició de la variabilitat
Mantenim la idea de la descomposició de la variabilitat total de les dades. En aquest cas la variabilitat total es descompon en: una part que representa la variabilitat deguda al factor A, una part que correspon a la variabilitat deguda al factor B, una part que correspon a la variabilitat deguda a la interacció dels factors A i B i finalment una part que és la variabilitat dintre dels tractaments o variabilitat del error o residual.
Posteriorment decidirem si els possibles efectes són significatius comparant les variabilitats dels efectes amb la variabilitat residual.
\[ \begin{aligned} SS_T = \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{ijk}-\overline{Y}_{...} )^2 = \\ bn \sum_{i=1}^{a} (\overline{Y}_{i..}-\overline{Y}_{...} )^2 + an \sum_{j=1}^{b} (\overline{Y}_{.j.}-\overline{Y}_{...} )^2 + \\ n \sum_{i=1}^{a} \sum_{j=1}^{b} (\overline{Y}_{ij.}-\overline{Y}_{i..}-\overline{Y}_{.j.}+\overline{Y}_{...} )^2 + \sum_{i=1}^{a} \sum_{j=1}^{b} \sum_{k=1}^{n} (Y_{ijk}-\overline{Y}_{ij.} )^2 = \\ SS_A + SS_B + SS_{AB} + SS_E \\ \end{aligned} \]
on
\(SS_T\) recull la variabilitat total de les dades,
\(SS_A\) recull la variabilitat deguda al factor A,
\(SS_B\) recull la variabilitat deguda al factor B,
\(SS_{AB}\) recull la variabilitat deguda a la interacció,
\(SS_E\) recull la variabilitat residual.
1.1.7 Esperances dels quadrats mitjans
1.1.7.1 Quadrats mitjans del error
Si anomenem \(MS_E = \frac{SS_E}{ab(n-1)}\) (quadrats mitjans de l’error) pot demostrar-se que sempre \(E(MS_E) = \sigma^2\).
1.1.7.2 Quadrats mitjans deguts al tractament A
Per una altra banda si anomenem \(MS_A = \frac{SS_A}{a-1}\) (quadrats mitjans del tractament A) pot demostrar-se que \(E(MS_A) = \sigma^2 + \frac{ bn \sum_{i=1}^a \alpha_i^2}{a-1}\).
Per tant i si la hipòtesi nul·la referent al factor \(A\) és certa \((H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0) \Rightarrow E(MS_A)=E(MS_E)=\sigma^2\) i per tant si plantegem
\[ F = \frac{MS_A}{MS_E} \approx 1. \]
En canvi si \(H_0\) és falsa la tendència ha de ser que \(F > 1\).
\[ F = \frac{MS_A}{MS_E} = \frac{\frac{SS_A}{a-1}}{\frac{SS_E}{ab(n-1)}}. \]
Si \(H_0\) és certa , \(F\) tendirà a ser propera a 1.
Si \(H_0\) és falsa , \(F\) tendirà a ser més gran que 1.
Per tant un criteri raonable per decidir consisteix en rebutjar \(H_0\) si \(F\) és prou gran.
Si \(F \ge c_\alpha \Rightarrow\) rebutgem \(H_0\), on \(c_\alpha\) és triat de forma que \(P_{H_0} (F \ge c_\alpha) = \alpha\).
Per trobar \(c_\alpha\) pot demostrar-se que sota la hipòtesi de errors normals i homocedàstics
\[ \epsilon_{ijk} \sim N(0,\sigma).\]
Sota la hipòtesi nul·la \(MS_A\) i \(MS_E\) segueixen una distribució \(\chi^2\) amb \(a-1\) i \(ab(n-1)\) graus de llibertat respectivament i són independents, i com a conseqüència
\[ F = \frac{MS_A}{MS_E} = \frac{\frac{SS_A}{a-1}}{\frac{SS_E}{ab(n-1)}} \sim F(a-1,ab(n-1))\]
on \(F(a-1,ab(n-1))\) és una distribució \(F\) de Fisher amb \(a-1\) graus de llibertat en el numerador i \(ab(n-1)\) graus de llibertat en el denominador.
1.1.7.3 Quadrats mitjans deguts al tractament B
Anàlogament si anomenem \(MS_B = \frac{SS_B}{b-1}\) pot demostrar-se que \(E(MS_B) = \sigma^2 + \frac{ an \sum_{j=1}^b \beta_j^2}{b-1}\).
Per tant i si la hipòtesi nul·la referent al factor \(B\) és certa \((H_0: \beta_1 = \beta_2 = \cdots = \beta_b = 0) \Rightarrow E(MS_B)=E(MS_E)=\sigma^2\) i per tant si plantegem
\[ F = \frac{MS_B}{MS_E} \approx 1 .\]
En canvi si \(H_0\) és falsa la tendència ha de ser que \(F > 1\).
Sota la hipòtesi nul·la \[ F = \frac{MS_B}{MS_E} = \frac{\frac{SS_B}{b-1}}{\frac{SS_E}{ab(n-1)}} \sim F(b-1,ab(n-1)).\]
1.1.7.4 Quadrats mitjans deguts a la interacció
Si anomenem \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) pot demostrar-se que \(E(MS_{AB}) = \sigma^2 + \frac{ n \sum_{i=1}^a \sum_{j=1}^b \gamma_{ij}^2}{(a-1)(b-1)}\).
Per tant i si la hipòtesi nul·la referent a la interacció és certa \(H_0: \gamma_{11} = \gamma_{12} = \cdots = \gamma_{ab} = 0 \Rightarrow E(MS_{AB})=E(MS_E)=\sigma^2\) i per tant si plantegem
\[ F = \frac{MS_{AB}}{MS_E} \approx 1. \]
En canvi si \(H_0\) és falsa la tendència ha de ser que \(F > 1\).
Sota la hipòtesi nul·la \[ F = \frac{MS_{AB}}{MS_E} = \frac{\frac{SS_{AB}}{(a-1)(b-1)}}{\frac{SS_E}{ab(n-1)}} \sim F((a-1)(b-1),ab(n-1)).\]
1.1.8 Taula de l’Anàlisi de la variància de dos factors (ANOVA)
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_E}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_E}\) |
| Interacció AB | \(SS_{AB}\) | \((a-1)(b-1)\) | \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) | \(F = \frac{MS_{AB}}{MS_E}\) |
| Error | \(SS_E\) | \(ab(n-1)\) | \(MS_E = \frac{SS_E}{ab(n-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
- Criteri de decisió
– Calcular el valor dels estadístics \(F\) de la taula a partir de les dades.
– Determinar els valors crítics \(F_{\alpha} (gl_{num},gl_{den})\) per a les corresponents distribucions \(F\) de Fisher d’acord amb \(\alpha\).
– Rebutjarem cada hipòtesi nul·la si l’estadístic \(F\) respectiu és més gran que el valor crític \(F \ge F_{\alpha} (gl_{num},gl_{den})\) corresponent o, equivalentment si obtenim el resultat amb un software estadístic, rebutjarem si el p-valor és menor que el nivell de significació.
1.1.9 Resolució amb R
La sintaxi vacuna*dosis representa la generació dels factors principals i la interacció.
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| vacuna | 3 | 140.557292 | 46.8524306 | 11.3339325 | 0.0000221 |
| dosis | 2 | 7.937917 | 3.9689583 | 0.9601189 | 0.3924335 |
| vacuna:dosis | 6 | 4.362083 | 0.7270139 | 0.1758698 | 0.9816806 |
| Residuals | 36 | 148.817500 | 4.1338194 | NA | NA |
En aquest cas direm que hi han diferències degudes a les vacunes, no a les dosis i tampoc hi ha una interacció significativa.
En forma alternativa, sense utilitzar l’asterisc i especificant tots els termes:
Gràfica descriptiva: Interaction plot
interaction.plot(immuno$vacuna,immuno$dosis,immuno$anticos,col='red',trace.label='dosis',xlab='Vacuna',ylab='Anticos')
1.1.9.1 Taules de coeficients i mitjanes
## Tables of means
## Grand mean
##
## 4.960417
##
## vacuna
## vacuna
## 1 2 3 4
## 7.542 5.150 4.342 2.808
##
## dosis
## dosis
## 1 2 3
## 5.000 4.444 5.437
##
## vacuna:dosis
## dosis
## vacuna 1 2 3
## 1 7.725 7.275 7.625
## 2 5.000 4.475 5.975
## 3 4.525 3.325 5.175
## 4 2.750 2.700 2.975## Tables of effects
##
## vacuna
## vacuna
## 1 2 3 4
## 2.5812 0.1896 -0.6187 -2.1521
##
## dosis
## dosis
## 1 2 3
## 0.0396 -0.5167 0.4771
##
## vacuna:dosis
## dosis
## vacuna 1 2 3
## 1 0.1438 0.2500 -0.3938
## 2 -0.1896 -0.1583 0.3479
## 3 0.1437 -0.5000 0.3563
## 4 -0.0979 0.4083 -0.3104- Algunes estimacions dels paràmetres del model:
– \(\hat{\mu} =\) 4.9604167
– \(\hat{\alpha}_1\) = 7.5416667 - 4.9604167= 2.58125
– \(\hat{\beta}_2\) = 4.44375 - 4.9604167= -0.5166667
– \(\hat{\gamma}_{23}\) = 5.975 - 5.15 - 5.4375 + 4.9604167= 0.3479167
– \(\hat{\sigma}^2\) = 4.1338194 \(\Rightarrow \hat{\sigma}\) = 2.0331796
1.1.10 Comparacions múltiples
Assumint que algun dels factors principals ha estat significatiu i la interacció no, podem plantejar-nos les comparacions múltiples per a detectar quins nivells del factor principal són diferents entre ells.
##
## Study: immuno.aov ~ "vacuna"
##
## HSD Test for anticos
##
## Mean Square Error: 4.133819
##
## vacuna, means
##
## anticos std r se Min Max Q25 Q50 Q75
## 1 7.541667 1.5382743 12 0.5869284 6.1 11.0 6.450 6.85 8.350
## 2 5.150000 1.8705857 12 0.5869284 2.9 10.2 4.125 5.15 5.600
## 3 4.341667 2.8833246 12 0.5869284 2.2 12.4 2.825 3.45 4.225
## 4 2.808333 0.6841828 12 0.5869284 1.8 3.8 2.275 2.70 3.375
##
## Alpha: 0.05 ; DF Error: 36
## Critical Value of Studentized Range: 3.808798
##
## Minimun Significant Difference: 2.235492
##
## Treatments with the same letter are not significantly different.
##
## anticos groups
## 1 7.541667 a
## 2 5.150000 b
## 3 4.341667 bc
## 4 2.808333 cEn aquest cas no seria correcte fer les comparacions múltiples pel factor dosis perquè no ha estat significatiu.
Però, com a exercici didàctic, què passa si ho fem:
##
## Study: immuno.aov ~ "dosis"
##
## HSD Test for anticos
##
## Mean Square Error: 4.133819
##
## dosis, means
##
## anticos std r se Min Max Q25 Q50 Q75
## 1 5.00000 2.196968 16 0.5082949 2.1 9.5 3.450 4.35 6.275
## 2 4.44375 2.027468 16 0.5082949 1.8 8.2 2.900 3.75 5.775
## 3 5.43750 3.262693 16 0.5082949 2.2 12.4 2.975 4.50 6.600
##
## Alpha: 0.05 ; DF Error: 36
## Critical Value of Studentized Range: 3.456758
##
## Minimun Significant Difference: 1.757053
##
## Treatments with the same letter are not significantly different.
##
## anticos groups
## 3 5.43750 a
## 1 5.00000 a
## 2 4.44375 aTots els resultats ens permeten reinterpretar el gràfic d’interaccions.
1.1.11 Interpretació de la interacció
A l’exemple anterior hem vist que la interacció no era significativa, això es traduïa en un paral·lelisme en el gràfic d’interaccions
Al següent apartat veurem un altre exemple força diferent.
1.1.11.1 Exemple 2: Grau de somnolència en fàrmacs analgèsics
En fàrmacs per reduir el dolor, el grau de somnolència és un dels efectes secundaris que de forma habitual s’avalua en base a la consideració del temps de reacció a un estímul (en segons). En aquesta direcció, estem interessats en contrastar l’efecte de tres medicaments (A, B i C) i en la forma d’administrar-los (C: càpsula, I: injectable). Els resultats obtinguts en individus de característiques similars són els següents
temps<-c(4.35,3.86,4.62,5.14,4.8,3.71,3.9,3.4,5.97,6.05,6.32,
6.07,1.35,1.62,0.34,0.63,1.08,1.51,1.61,0.94,0.84,0.81,1.67,0.96)
admin<-gl(2,12)
farmac<-factor(rep(rep(1:3,each=4),2))
somnol<-data.frame(temps,admin,farmac)
interaction.plot(somnol$farmac,somnol$admin,somnol$temps,col='red',trace.label = 'admin')
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| farmac | 2 | 4.449900 | 2.2249500 | 10.27420 | 0.0010553 |
| admin | 1 | 83.738704 | 83.7387042 | 386.68215 | 0.0000000 |
| farmac:admin | 2 | 5.749633 | 2.8748167 | 13.27511 | 0.0002869 |
| Residuals | 18 | 3.898025 | 0.2165569 | NA | NA |
Quan la interacció és significativa això ens vol dir que els efectes principals no són independents un de l’altre i per tant es fa difícil d’interpretar el sentit de les possibles diferències entre els nivells dels factors. Es més recomanable analitzar les diferències entre els nivells de la interacció o fer comparacions per a un factor dins cada nivell de l’altre factor.
##
## Study: somnol.aov ~ c("admin", "farmac")
##
## HSD Test for temps
##
## Mean Square Error: 0.2165569
##
## admin:farmac, means
##
## temps std r se Min Max Q25 Q50 Q75
## 1:1 4.4925 0.5341270 4 0.2326784 3.86 5.14 4.2275 4.485 4.7500
## 1:2 3.9525 0.6014081 4 0.2326784 3.40 4.80 3.6325 3.805 4.1250
## 1:3 6.1025 0.1512999 4 0.2326784 5.97 6.32 6.0300 6.060 6.1325
## 2:1 0.9850 0.5995832 4 0.2326784 0.34 1.62 0.5575 0.990 1.4175
## 2:2 1.2850 0.3252179 4 0.2326784 0.94 1.61 1.0450 1.295 1.5350
## 2:3 1.0700 0.4052160 4 0.2326784 0.81 1.67 0.8325 0.900 1.1375
##
## Alpha: 0.05 ; DF Error: 18
## Critical Value of Studentized Range: 4.49442
##
## Minimun Significant Difference: 1.045754
##
## Treatments with the same letter are not significantly different.
##
## temps groups
## 1:3 6.1025 a
## 1:1 4.4925 b
## 1:2 3.9525 b
## 2:2 1.2850 c
## 2:3 1.0700 c
## 2:1 0.9850 cPodem veure diverses situacions a la següent figura extreta del llibre de G. Quinn i M. Keough Experimental Design and Data Analysis for Biologists, Cambridge University Press, 2002.
1.1.12 Verificació de les premisses del ANOVA
1.1.12.1 Exemple 2: Grau de somnolència en fàrmacs analgèsics
1.1.12.1.1 Normalitat dels residus
##
## Shapiro-Wilk normality test
##
## data: somnol.aov$residuals
## W = 0.95564, p-value = 0.3572
1.1.12.2 Exemple 1: vacunes i dosis
1.1.12.2.1 Normalitat dels residus
##
## Shapiro-Wilk normality test
##
## data: immuno.aov$residuals
## W = 0.87414, p-value = 0.0001017
1.1.12.2.2 Homogeneïtat de variàncies
##
## Bartlett test of homogeneity of variances
##
## data: anticos by interaction(vacuna, dosis)
## Bartlett's K-squared = 26.217, df = 11, p-value = 0.006025
En aquest cas veiem clarament tant una desviació de la normalitat com una heteroscedasticitat (variàncies desiguals). Una possible forma de corregir el problema pot ser una transformació de la variable resposta.
Algunes de les transformacions més usuals poden trobar-se a la taula següent:
| Tipus de dades | Transformació | Instrucció R |
|---|---|---|
| Mesures continues | ln(x) | log(x) |
| - | ln(x+1) | log(x+1) |
| - | log(x) | log(x,10) |
| Comptatges | \(\sqrt{x}\) | sqrt(x) |
| Percentatges | arcsin | asin(sqrt(x))*180/pi |
1.1.12.2.3 Anàlisi amb les dades transformades
Nosaltres hem optat per aplicar la transformació logarítmica en base \(e\).
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| vacuna | 3 | 6.4060309 | 2.1353436 | 16.3520779 | 0.0000007 |
| dosis | 2 | 0.2059845 | 0.1029923 | 0.7886962 | 0.4621314 |
| vacuna:dosis | 6 | 0.1281582 | 0.0213597 | 0.1635688 | 0.9848036 |
| Residuals | 36 | 4.7010766 | 0.1305855 | NA | NA |
##
## Shapiro-Wilk normality test
##
## data: immuno.aov.trans$residuals
## W = 0.94396, p-value = 0.0231
##
## Bartlett test of homogeneity of variances
##
## data: log(anticos) by interaction(vacuna, dosis)
## Bartlett's K-squared = 15.893, df = 11, p-value = 0.1451
Com veiem els problemes de no normalitat i hetersocedasticitat s’han corregit.
1.1.13 Cas d’una rèplica per nivell de la interacció
Si tenim una única rèplica per cada combinació de nivells dels dos factors (una rèplica per casella) \(n=1\) la suma de quadrats del error \(SS_E\) té 0 graus de llibertat i \(\sigma^2\) no és estimable.
La solució es no incloure la interacció en el model. L’equació del model passa a ser:
\[ Y_{ij} = \mu + \alpha_i + \beta_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b , \]
i en aquest cas \(SS_T = SS_A + SS_B + SS_E\), adaptant les fórmules al cas \(n=1\).
Formalment la interacció no pot separar-se dels residus.
La taula ANOVA queda
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_E}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_E}\) |
| Error | \(SS_E\) | \((a-1)(b-1)\) | \(MS_E = \frac{SS_E}{(a-1)(b-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
Per verificar en el cas d’una rèplica si és assumible la suposició de no interacció existeixen alguns tests com per exemple el Test de additivitat de Tukey.
1.1.13.1 Subexemple de l’exemple 1: vacunes i dosis
anticos<-c(6.1,8.2,6.3,4.3,5.2,5.6,4.0,4.9,3.1,2.1,3.7,2.5)
vacuna<-gl(4,3)
dosis<-factor(rep(1:3,times=4))
immuno1<-data.frame(anticos,vacuna,dosis)
## anova amb la interacció
immuno1.aov<-aov(anticos~vacuna*dosis,immuno1)
anova(immuno1.aov)| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| vacuna | 3 | 27.086667 | 9.0288889 | NaN | NaN |
| dosis | 2 | 4.291667 | 2.1458333 | NaN | NaN |
| vacuna:dosis | 6 | 2.288333 | 0.3813889 | NaN | NaN |
| Residuals | 0 | 0.000000 | NaN | NA | NA |
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| vacuna | 3 | 27.086667 | 9.0288889 | 23.673707 | 0.0010034 |
| dosis | 2 | 4.291667 | 2.1458333 | 5.626366 | 0.0420611 |
| Residuals | 6 | 2.288333 | 0.3813889 | NA | NA |
1.2 Dos factors aleatoris
1.2.1 Exemple 3. Variabilitat de la mostra i del analista.
Es vol avaluar la variabilitat en la determinació de la qualitat d’un medi de cultiu per tal d’estudiar la rapidesa d’aparició de noves colònies d’Escherichia coli, en base al recompte d’individus en un cultiu durant un període de 48 hores. Es creu que els errors en el comptatge es poden atribuir a la preparació de la mostra o al comptatge per part dels analistes. Per aquest motiu, s’han realitzat 3 preparacions triades a l’atzar de la mateixa base i es realitza l’anàlisi de cada una d’elles per part de dos operaris triats també a l’atzar. Volem valorar la variabilitat aportada a la variable resposta pels diferents factors.
Els resultats obtinguts són:
compt<-c(343,344,341,339,365,364,359,361,321,325,317,322)
prep<-gl(3,4)
analist<-factor(rep(rep(1:2,each=2),3))
dades<-data.frame(prep,analist,compt)
dades| prep | analist | compt |
|---|---|---|
| 1 | 1 | 343 |
| 1 | 1 | 344 |
| 1 | 2 | 341 |
| 1 | 2 | 339 |
| 2 | 1 | 365 |
| 2 | 1 | 364 |
| 2 | 2 | 359 |
| 2 | 2 | 361 |
| 3 | 1 | 321 |
| 3 | 1 | 325 |
| 3 | 2 | 317 |
| 3 | 2 | 322 |
L’objectiu es veure si l’analista o la preparació de la mostra tenen influència en la variable resposta.
No estem interessants especialment ni en els dos analistes en concret que han fet el recompte ni en les tres preparacions analitzades, voldríem que les conclusions fossin vàlides pel conjunt de preparacions i analistes possibles, els dos factors són aleatoris.
1.2.2 Model lineal del ANOVA de dos factors aleatoris
Suposem ara que tenim un experiment amb dues fonts de variabilitat controlades \(A\) i \(B\) però que les dues són aleatòries, es a dir els seus nivells són mostres aleatòries de mida \(a\) i \(b\) respectivament d’un conjunt més gran de possibles nivells.
El model lineal és ara
\[ Y_{ijk} = \mu + A_i + B_j + AB_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]
on
\(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\).
\(\mu\) és la mitjana general.
\(A_i\) és una variable aleatòria, no un paràmetre constant. Assumim que \(A_i \sim N(0,\sigma_A)\).
\(B_j\) és una variable aleatòria, no un paràmetre constant. Assumim que \(B_j \sim N(0,\sigma_B)\).
\(AB_{ij}\) és una variable aleatòria, no un paràmetre constant. Assumim que \(AB_{ij} \sim N(0,\sigma_{AB})\).
\(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ijk} \sim N(0,\sigma)\).
Assumirem també que les v.a. \(A_i\), \(B_j\), \(AB_{ij}\) i \(\epsilon_{ijk}\) són v.a. independents.
i per tant cada observació experimental
\[ Y_{ijk} \sim N(\mu,\sqrt{\sigma^2_A+\sigma^2_B+\sigma^2_{AB}+\sigma^2} ).\]
1.2.2.1 Dependència entre les observacions
Tal i com passa en el model aleatori d’un factor, l’existència de factors aleatoris introdueix una dependència entre les observacions dintre d’algunes cel·les del disseny. Pot demostrar-se que:
\[ Cov(Y_{ijk},Y_{i'j'k'}) = \left\{ \begin{array}{ll} \sigma^2_A+\sigma^2_B+\sigma^2_{AB} & \mbox{ si } i=i', j=j', k \ne k' \\ \sigma^2_A & \mbox{ si } i=i', j \ne j' \\ \sigma^2_B & \mbox{ si } i \ne i', j = j' \\ 0 & \mbox{ si } i \ne i', j \ne j' \end{array} \right. \]
Evidentment si \(i=i', j=j', k = k'\) llavors
\[ Cov(Y_{ijk},Y_{ijk})=Var(Y_{ijk})=\sigma^2_A+\sigma^2_B+\sigma^2_{AB}+\sigma^2\]
1.2.3 Contrasts d’hipòtesis
Els contrasts de més interès ara són els següents
Variabilitat aportada pel factor A \[ H_0: \sigma_A^2 = 0 \] \[ H_1: \sigma_A^2 > 0\]
Variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]
Variabilitat aportada per la interacció \[ H_0: \sigma_{AB}^2 = 0 \] \[ H_1: \sigma_{AB}^2 > 0\]
La descomposició de la suma de quadrats continua sent vàlida \(SS_A + SS_B +SS_{AB} + SS_E\) i les sumes de quadrats es calculen de la mateixa forma, però els quadrats mitjans tenen un altre sentit i les esperances són ara les següents
- \(E(MS_E) = \sigma^2\)
- \(E(MS_A) = \sigma^2 + bn \sigma^2_A+ n \sigma^2_{AB}\)
- \(E(MS_B) = \sigma^2 + an \sigma^2_B+ n \sigma^2_{AB}\)
- \(E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB}\)
Per tant els estadístics F per decidir en cadascú dels contrasts anteriors són
\[ F = \frac{MS_A}{MS_{AB}} \sim F(a-1,(a-1)(b-1)) \]
\[ F = \frac{MS_B}{MS_{AB}} \sim F(b-1,(a-1)(b-1)) \]
\[ F = \frac{MS_{AB}}{MS_E} \sim F((a-1)(b-1),ab(n-1)) \]
i la taula ANOVA és la següent
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_{AB}}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_{AB}}\) |
| Interacció AB | \(SS_{AB}\) | \((a-1)(b-1)\) | \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) | \(F = \frac{MS_{AB}}{MS_E}\) |
| Error | \(SS_E\) | \(ab(n-1)\) | \(MS_E = \frac{SS_E}{ab(n-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
1.2.4 Resolució de l’exemple 3
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| prep | 2 | 3362.0000000 | 1681.0000000 | 395.5294118 | 0.0000004 |
| analist | 1 | 44.0833333 | 44.0833333 | 10.3725490 | 0.0181225 |
| prep:analist | 2 | 0.6666667 | 0.3333333 | 0.0784314 | 0.9254977 |
| Residuals | 6 | 25.5000000 | 4.2500000 | NA | NA |
Els resultats de l’ANOVA encara no són correctes. R assumeix per defecte que tots els factors són fixos.
1.2.4.1 Modificació manual de la taula ANOVA
## Modificació de la taula a ma
taula<-anova(res.aov)
taula[1,4]<-taula[1,3]/taula[3,3]
taula[2,4]<-taula[2,3]/taula[3,3]
taula[1,5]<-1-pf(taula[1,4],taula[1,1],taula[3,1])
taula[2,5]<-1-pf(taula[2,4],taula[2,1],taula[3,1])
taula| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| prep | 2 | 3362.0000000 | 1681.0000000 | 5.04300e+03 | 0.0001983 |
| analist | 1 | 44.0833333 | 44.0833333 | 1.32250e+02 | 0.0074767 |
| prep:analist | 2 | 0.6666667 | 0.3333333 | 7.84314e-02 | 0.9254977 |
| Residuals | 6 | 25.5000000 | 4.2500000 | NA | NA |
1.2.5 Estimació dels paràmetres del model. Components de la variància.
\(\hat{\mu} = \overline{Y}_{...}\)
\(\hat{\sigma}^2 = MS_E\)
\(\hat{\sigma}^2_A = \frac{MS_A-MS_{AB}}{bn}\)
\(\hat{\sigma}^2_B = \frac{MS_B-MS_{AB}}{an}\)
\(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n}\)
Aquestes estimacions s’obtenen pel mètode dels moments a partir de les esperances dels quadrats mitjans
\[ E(MS_E) = \sigma^2 \,\,\,\, ; \,\,\,\, E(MS_A) = \sigma^2 + bn \sigma^2_A+ n \sigma^2_{AB} \] \[ E(MS_B) = \sigma^2 + an \sigma^2_B+ n \sigma^2_{AB} \,\,\,\, ; \,\,\,\, E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB} \]
1.2.6 Estimació de components de la variància a l’exemple 3
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| prep | 2 | 3362.0000000 | 1681.0000000 | 5.04300e+03 | 0.0001983 |
| analist | 1 | 44.0833333 | 44.0833333 | 1.32250e+02 | 0.0074767 |
| prep:analist | 2 | 0.6666667 | 0.3333333 | 7.84314e-02 | 0.9254977 |
| Residuals | 6 | 25.5000000 | 4.2500000 | NA | NA |
\(\hat{\sigma}^2 = MS_E =\) 4.25.
\(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n} =\) -1.9583333 \(\Rightarrow 0\).
\(\hat{\sigma}^2_A = \frac{MS_A-MS_{AB}}{bn} =\) 420.1666667.
\(\hat{\sigma}^2_B = \frac{MS_B-MS_{AB}}{an} =\) 7.2916667.
compA<-(taula[1,3]-taula[3,3])/(4)
compB<-(taula[2,3]-taula[3,3])/(6)
compAB<-(taula[3,3]-taula[4,3])/(2)
compE<-taula[4,3]
if (compA<0){compA<-0}
if (compAB<0){compAB<-0}
if (compB<0){compB<-0}
perA<-100*compA/(compA+compB+compAB+compE)
perB<-100*compB/(compA+compB+compAB+compE)
perAB<-100*compAB/(compA+compB+compAB+compE)
perE<-100*compE/(compA+compB+compAB+compE)
c(compA,compB,compAB,compE)## [1] 420.166667 7.291667 0.000000 4.250000## [1] 97.326513 1.689026 0.000000 0.9844611.2.6.1 Resolució amb el package mixlm
## Analysis of variance (unrestricted model)
## Response: compt
## Mean Sq Sum Sq Df F value Pr(>F)
## prep 1681.00 3362.00 2 5043.00 0.0002
## analist 44.08 44.08 1 132.25 0.0075
## prep:analist 0.33 0.67 2 0.08 0.9255
## Residuals 4.25 25.50 6 - -
##
## Err.term(s) Err.df VC(SS)
## 1 prep (3) 2 420.17
## 2 analist (3) 2 7.29
## 3 prep:analist (4) 6 -1.96
## 4 Residuals - - 4.25
## (VC = variance component)
##
## Expected mean squares
## prep (4) + 2 (3) + 4 (1)
## analist (4) + 2 (3) + 6 (2)
## prep:analist (4) + 2 (3)
## Residuals (4)1.2.7 Resolució alternativa de l’exemple 3 amb el package lme4
Ho farem només amb el package lme4, nlme amb la funció lme() sempre assumeix que els factors estan jerarquitzats (niats) i és molt difícil modelar factors aleatoris creuats.
1.2.7.1 Package lme4
library(lme4)
library(lmerTest)
res.lmer<-lmer(compt~1+(1|prep)+(1|analist)+(1|prep:analist))
summary(res.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: compt ~ 1 + (1 | prep) + (1 | analist) + (1 | prep:analist)
##
## REML criterion at convergence: 61.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3909 -0.6030 0.1357 0.5627 1.3738
##
## Random effects:
## Groups Name Variance Std.Dev.
## prep:analist (Intercept) 0.000 0.000
## prep (Intercept) 419.435 20.480
## analist (Intercept) 6.802 2.608
## Residual 3.271 1.809
## Number of obs: 12, groups: prep:analist, 6; prep, 3; analist, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 341.750 11.979 2.095 28.53 0.000949 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')## 2.5 % 97.5 %
## .sig01 0.0000000 3.426253
## .sig02 8.8870627 49.709260
## .sig03 0.6087447 18.734220
## .sigma 1.1866227 3.247714
## (Intercept) 314.4782199 369.021888| npar | logLik | AIC | LRT | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|
| 5 | -30.91114 | 71.82228 | NA | NA | NA | |
| (1 | prep) | 4 | -37.50172 | 83.00345 | 13.181166 | 1 | 0.0002828 |
| (1 | analist) | 4 | -32.67240 | 73.34479 | 3.522511 | 1 | 0.0605407 |
| (1 | prep:analist) | 4 | -30.91114 | 69.82228 | 0.000000 | 1 | 1.0000000 |
1.2.8 Dos factors aleatoris amb una observació per casella
En aquest cas no podem incloure la interacció en el model i el model que resta és
\[ Y_{ij} = \mu + A_i + B_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \] on
\(\mu\) és la mitjana general,
\(A_i \sim N(0,\sigma_A)\),
\(B_j \sim N(0,\sigma_B)\),
\(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ij} \sim N(0,\sigma)\).
Assumirem també que les v.a. \(A_i\), \(B_j\) i \(\epsilon_{ijk}\) són v.a. independents.
En aquest cas \(SS_T = SS_A + SS_B + SS_E\) adaptant les fórmules al cas \(n=1\).
Els únics contrasts són:
Variabilitat aportada pel factor A \[ H_0: \sigma_A^2 = 0 \] \[ H_1: \sigma_A^2 > 0\]
Variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]
La taula ANOVA queda
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_E}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_E}\) |
| Error | \(SS_E\) | \((a-1)(b-1)\) | \(MS_E = \frac{SS_E}{(a-1)(b-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
1.2.8.1 Estimacions dels paràmetres del model
\(\hat{\mu} = \overline{Y}_{...}\)
\(\hat{\sigma}^2 = MS_E\)
\(\hat{\sigma}^2_A = \frac{MS_A-MS_E}{b}\)
\(\hat{\sigma}^2_B = \frac{MS_B-MS_E}{a}\)
1.3 Un factor fix i un factor aleatori
1.3.1 Exemple 4. Estudi d’adhesió cel·lular
En un estudi d’adhesió cel·lular es van analitzar 4 tipus cel·lulars segons la seva capacitat d’adhesió a un determinat substrat. Per raons logístiques es van realitzar tres experiments independents en tres dies diferents. En cada experiment es van analitzar els 4 tipus cel·lulars amb un total de 5 rèpliques per cada tipus. Volem valorar la influència en la variable resposta dels diferents factors que intervenen a l’estudi.
Els resultats són els següents
adhesio1<-c(0.89,0.95,1.31,1.01,1.03,1.24,0.99,0.65,0.76,0.81,
0.94,0.77,0.78,0.72,0.59,1.36,1.1,0.61,1.14,0.91)
adhesio2<-c(0.69,0.86,0.69,0.75,0.77,0.46,0.97,0.83,0.57,1.05,
0.69,0.54,0.78,0.78,0.76,1.04,1.01,0.89,0.84,0.82)
adhesio3<-c(0.36,0.46,0.25,0.93,0.74,0.18,0.32,0.31,0.56,0.38,
0.25,0.19,0.19,0.47,0.15,0.8,0.19,0.62,0.52,0.54)
adhesio<-c(adhesio1,adhesio2,adhesio3)
cell<-rep(gl(4,5),3)
experiment<- gl(3,20)
myData <- data.frame(experiment,cell,adhesio)El factor important és el tipus cel·lular, factor fix amb 4 nivells, però també volem estudiar la variabilitat que pot introduir el dia en que es fa l’experiment, per tant incluïm el factor experiment en el model. Tenim per tant un factor fix i un factor aleatori.
interaction.plot(myData$experiment,myData$cell,response=myData$adhesio,main='Interaction Plot',col='red',trace.label = 'cell',xlab='Experiment',ylab='Adhesió')
1.3.2 Dissenys mixtos. Un factor fix i un factor aleatori
Suposem ara que tenim un experiment amb dues fonts de variabilitat controlades \(A\) i \(B\) però on una es fixa \(A\) i l’altre aleatòria \(B\), es a dir els seus nivells són una mostra aleatòria de mida \(b\) d’un conjunt més gran de possibles nivells.
El model lineal és ara
\[ Y_{ijk} = \mu + \alpha_i + B_j + (\alpha B)_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]
on
\(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),
\(\mu\) és la mitjana general,
\(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A,
\(B_j\) és una variable aleatòria, no un paràmetre constant. Assumim que \(B_j \sim N(0,\sigma_B)\),
\((\alpha B)_{ij}\) és una variable aleatòria, no un paràmetre constant. Assumim que \((\alpha B)_{ij} \sim N(0,\sigma_{AB})\). (Tota interacció on al menys un dels termes sigui aleatori és aleatòria).
Assumim també la restricció \(\sum_{i=1}^a \alpha_i = 0\).
\(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ijk} \sim N(0,\sigma)\).
Assumirem també que les v.a. \(B_j\), \((\alpha B)_{ij}\) i \(\epsilon_{ijk}\) són v.a. independents.
i per tant cada observació experimental
\[ Y_{ijk} \sim N(\mu+\alpha_i,\sqrt{\sigma^2_B+\sigma^2_{AB}+\sigma^2} )\]
1.3.2.1 Dependència entre les observacions
\[ Cov(Y_{ijk},Y_{i'j'k'}) = \left\{ \begin{array}{ll} \sigma^2_B+\sigma^2_{AB} & \mbox{ si } i=i', j=j', k \ne k' \\ \sigma^2_B & \mbox{ si } i \ne i', j = j' \\ 0 & \mbox{ si } i \ne i', j \ne j' \end{array} \right. \]
Evidentment si \(i=i', j=j', k = k'\) llavors
\[ Cov(Y_{ijk},Y_{ijk})=Var(Y_{ijk})=\sigma^2_B+\sigma^2_{AB}+\sigma^2\]
1.3.3 Dissenys mixtos. Contrasts d’hipòtesis
Els contrasts de més interès ara són els següents
Hipòtesi sobre la significació del factor A \[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]
Hipòtesi sobre la variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]
Hipòtesi sobre la variabilitat aportada per la interacció \[ H_0: \sigma_{AB}^2 = 0 \] \[ H_1: \sigma_{AB}^2 > 0\]
Les esperances dels quadrats mitjans són les següents
- \(E(MS_E) = \sigma^2\)
- \(E(MS_A) = \sigma^2 + n \sigma^2_{AB}+ \frac{ bn \sum_{i=1}^a \alpha_i^2}{a-1}\)
- \(E(MS_B) = \sigma^2 + an \sigma^2_B+ n \sigma^2_{AB}\)
- \(E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB}\)
Per tant els estadístics F per decidir en cadascú dels contrasts anteriors són
\[ F = \frac{MS_A}{MS_{AB}} \sim F(a-1,(a-1)(b-1)) \]
\[ F = \frac{MS_B}{MS_{AB}} \sim F(b-1,(a-1)(b-1)) \]
\[ F = \frac{MS_{AB}}{MS_E} \sim F((a-1)(b-1),ab(n-1)) \]
i la taula ANOVA és la següent:
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_{AB}}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_{AB}}\) |
| Interacció AB | \(SS_{AB}\) | \((a-1)(b-1)\) | \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) | \(F = \frac{MS_{AB}}{MS_E}\) |
| Error | \(SS_E\) | \(ab(n-1)\) | \(MS_E = \frac{SS_E}{ab(n-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
1.3.3.1 Estimació dels paràmetres del model. Components de la variància.
\(\hat{\mu} = \overline{Y}_{...}\)
\(\hat{\alpha_i} = \overline{Y}_{i..} - \overline{Y}_{...} \,\,\,\, (\alpha_i = \mu_i - \mu)\)
\(\hat{\sigma}^2 = MS_E\)
\(\hat{\sigma}^2_B = \frac{MS_B-MS_{AB}}{an}\)
\(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n}\)
1.3.4 Dissenys mixtos. Model restringit
Suposem ara que tenim un experiment amb dues fonts de variabilitat controlades \(A\) i \(B\) però on una es fixa \(A\) i l’altre aleatòria \(B\), es a dir els seus nivells són una mostra aleatòria de mida \(b\) d’un conjunt més gran de possibles nivells.
El model lineal és
\[ Y_{ijk} = \mu + \alpha_i + B_j + (\alpha B)_{ij} +\epsilon_{ijk} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,n_{ij} \]
Però ara modifiquem les restriccions als paràmetres on
\(Y_{ijk}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),
\(\mu\) és una mitjana general,
\(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A amb \(\sum_{i=1}^a \alpha_i = 0\),
\(B_j \sim N(0,\sigma_B)\),
\(\epsilon_{ijk}\) és l’error aleatori corresponent a l’observació \(i,j,k\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ijk} \sim N(0,\sigma)\).
La nova restricció s’imposa sobre \((\alpha B)_{ij}\):
\[ (\alpha B)_{ij} \sim N \left( 0 , \sqrt{\frac{a-1}{a}}\sigma_{AB} \right) \,\,\,\, , \,\,\,\, \sum_{i=1}^a (\alpha B)_{ij} = 0, \,\,\, j=1,...,b \]
A causa d’aquesta nova restricció el model es coneix com model restringit. La forma de la variància de \((\alpha B)_{ij}\) simplifica algunes expressions de les esperances de les sumes de quadrats.
1.3.5 Model restringit. Contrasts d’hipòtesis
Hipòtesi sobre la significació del factor A \[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]
Hipòtesi sobre la variabilitat aportada pel factor B \[ H_0: \sigma_B^2 = 0 \] \[ H_1: \sigma_B^2 > 0\]
Hipòtesi sobre la variabilitat aportada per la interacció \[ H_0: \sigma_{AB}^2 = 0 \] \[ H_1: \sigma_{AB}^2 > 0\]
Les esperances dels quadrats mitjans són les següents
- \(E(MS_E) = \sigma^2\)
- \(E(MS_A) = \sigma^2 + n \sigma^2_{AB}+ \frac{ bn \sum_{i=1}^a \alpha_i^2}{a-1}\)
- \(E(MS_B) = \sigma^2 + an \sigma^2_B\)
- \(E(MS_{AB}) = \sigma^2 + n \sigma^2_{AB}\)
Per tant els estadístics F per decidir en cadascú dels contrasts anteriors són
\[ F = \frac{MS_A}{MS_{AB}} \sim F(a-1,(a-1)(b-1)) \]
\[ F = \frac{MS_B}{MS_{E}} \sim F(b-1,ab(n-1)) \]
\[ F = \frac{MS_{AB}}{MS_E} \sim F((a-1)(b-1),ab(n-1)) \]
i la taula ANOVA és la següent
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_{AB}}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_{E}}\) |
| Interacció AB | \(SS_{AB}\) | \((a-1)(b-1)\) | \(MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}\) | \(F = \frac{MS_{AB}}{MS_E}\) |
| Error | \(SS_E\) | \(ab(n-1)\) | \(MS_E = \frac{SS_E}{ab(n-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
1.3.5.1 Estimació dels paràmetres del model. Components de la variància.
\(\hat{\mu} = \overline{Y}_{...}\)
\(\hat{\alpha_i} = \overline{Y}_{i..} - \overline{Y}_{...} \,\,\,\, (\alpha_i = \mu_i - \mu)\)
\(\hat{\sigma}^2 = MS_E\)
\(\hat{\sigma}^2_B = \frac{MS_B-MS_{E}}{an}\)
\(\hat{\sigma}^2_{AB} = \frac{MS_{AB}-MS_E}{n}\)
1.3.6 Restringit vs. no restringit. Quin model utilitzar?
Algunes persones i alguns paquets de software tenen preferències per una opció per defecte i l’apliquen sempre.
Però ens podem preguntar si hi han situacions on un model s’ajusta millor que l’altre.
La conseqüència més immediata de la restricció imposada al model restringit és que
\[ Cov ((\alpha B)_{ij},(\alpha B)_{i'j}) = -\sigma^2_{AB}/a. \]
Es a dir els termes d’interacció dintre de cada nivell del factor aleatori estan correlacionats negativament.
Això es tradueix en la següent estructura de dependència entre observacions pertanyents al mateix nivell del factor aleatori en el model restringit:
\[ Cov (Y_{ijk},Y_{i'jk'}) = \sigma_B^2-\frac{\sigma_{AB}^2}{a}. \]
Mentre que en el model no restringit
\[ Cov (Y_{ijk},Y_{i'jk'}) = \sigma_B^2. \]
Com a detall important en el model restringit la covariància pot ser negativa o positiva mentre que en el model no restringit la covariància sempre és positiva.
- Una situació en que això es reflecteix és en experiments amb recursos limitats.
Per exemple imaginem que estem comparant el creixement de dues espècies de plantes (factor fix) i que la parcel·la on creixen és un factor aleatori de b nivells al tractar-se de parcel·les seleccionades aleatòriament. Cada parcel·la es divideix en 2n unitats experimentals i es produeix una assignació aleatòria dels individus de forma que cada espècie té n individus en cada parcel·la. Degut a que els recursos (nutrients, aigua, espai) a cada parcel·la son limitats, esperem correlacions negatives entre els pesos finals de les plantes de cada parcel·la. Per tant en aquest cas semblaria millor el model restringit.
1.3.7 Resolució de l’exemple 4
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| cell | 3 | 0.5753467 | 0.1917822 | 5.3827311 | 0.0028246 |
| experiment | 2 | 2.7526633 | 1.3763317 | 38.6293533 | 0.0000000 |
| cell:experiment | 6 | 0.1201633 | 0.0200272 | 0.5621019 | 0.7581959 |
| Residuals | 48 | 1.7102000 | 0.0356292 | NA | NA |
1.3.7.1 Resolució manual del model no restringit
taula<-anova(mixed.aov)
taula[1,4]<-taula[1,3]/taula[3,3]
taula[2,4]<-taula[2,3]/taula[3,3]
taula[1,5]<-1-pf(taula[1,4],taula[1,1],taula[3,1])
taula[2,5]<-1-pf(taula[2,4],taula[2,1],taula[3,1])
taula| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| cell | 3 | 0.5753467 | 0.1917822 | 9.5760770 | 0.0105230 |
| experiment | 2 | 2.7526633 | 1.3763317 | 68.7230436 | 0.0000732 |
| cell:experiment | 6 | 0.1201633 | 0.0200272 | 0.5621019 | 0.7581959 |
| Residuals | 48 | 1.7102000 | 0.0356292 | NA | NA |
1.3.7.2 Resolució amb el package mixlm
## Analysis of variance (unrestricted model)
## Response: adhesio
## Mean Sq Sum Sq Df F value Pr(>F)
## cell 0.19 0.58 3 9.58 0.0105
## experiment 1.38 2.75 2 68.72 0.0001
## cell:experiment 0.02 0.12 6 0.56 0.7582
## Residuals 0.04 1.71 48 - -
##
## Err.term(s) Err.df VC(SS)
## 1 cell (3) 6 fixed
## 2 experiment (3) 6 0.06782
## 3 cell:experiment (4) 48 -0.00312
## 4 Residuals - - 0.03563
## (VC = variance component)
##
## Expected mean squares
## cell (4) + 5 (3) + 15 Q[1]
## experiment (4) + 5 (3) + 20 (2)
## cell:experiment (4) + 5 (3)
## Residuals (4)1.3.7.3 Resolució manual del model restringit
taula<-anova(mixed.aov)
taula[1,4]<-taula[1,3]/taula[3,3]
taula[1,5]<-1-pf(taula[1,4],taula[1,1],taula[3,1])
taula| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| cell | 3 | 0.5753467 | 0.1917822 | 9.5760770 | 0.0105230 |
| experiment | 2 | 2.7526633 | 1.3763317 | 38.6293533 | 0.0000000 |
| cell:experiment | 6 | 0.1201633 | 0.0200272 | 0.5621019 | 0.7581959 |
| Residuals | 48 | 1.7102000 | 0.0356292 | NA | NA |
1.3.7.4 Resolució amb el package mixlm
## Analysis of variance (restricted model)
## Response: adhesio
## Mean Sq Sum Sq Df F value Pr(>F)
## cell 0.19 0.58 3 9.58 0.0105
## experiment 1.38 2.75 2 38.63 0.0000
## cell:experiment 0.02 0.12 6 0.56 0.7582
## Residuals 0.04 1.71 48 - -
##
## Err.term(s) Err.df VC(SS)
## 1 cell (3) 6 fixed
## 2 experiment (4) 48 0.06704
## 3 cell:experiment (4) 48 -0.00312
## 4 Residuals - - 0.03563
## (VC = variance component)
##
## Expected mean squares
## cell (4) + 5 (3) + 15 Q[1]
## experiment (4) + 20 (2)
## cell:experiment (4) + 5 (3)
## Residuals (4)1.3.8 Comparacions múltiples
L’única part complicada dels contrasts o comparacions per parelles sobre els factors principals fixos és assegurar-se que s’utilitza el terme d’error correcte si el model conté factors aleatoris. En aquestes situacions, el terme d’error per al factor fix, i per tant per qualsevol contrast sobre el valor mitjà d’aquests factors, sol ser la interacció en lloc del residual.
df<-taula[3,1]
MSerror<-taula[3,3]
HSD.test(myData$adhesio,myData$cell,df,MSerror,console=T,group=F)##
## Study: myData$adhesio ~ myData$cell
##
## HSD Test for myData$adhesio
##
## Mean Square Error: 0.02002722
##
## myData$cell, means
##
## myData.adhesio std r se Min Max Q25 Q50 Q75
## 1 0.7793333 0.2731684 15 0.03653968 0.25 1.31 0.690 0.77 0.940
## 2 0.6720000 0.3114299 15 0.03653968 0.18 1.24 0.420 0.65 0.900
## 3 0.5733333 0.2619887 15 0.03653968 0.15 0.94 0.360 0.69 0.775
## 4 0.8260000 0.2951465 15 0.03653968 0.19 1.36 0.615 0.84 1.025
##
## Alpha: 0.05 ; DF Error: 6
## Critical Value of Studentized Range: 4.895599
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## 1 - 2 0.10733333 0.2605 -0.07155029 0.28621696
## 1 - 3 0.20600000 0.0277 * 0.02711638 0.38488362
## 1 - 4 -0.04666667 0.8043 -0.22555029 0.13221696
## 2 - 3 0.09866667 0.3163 -0.08021696 0.27755029
## 2 - 4 -0.15400000 0.0881 . -0.33288362 0.02488362
## 3 - 4 -0.25266667 0.0109 * -0.43155029 -0.073783041.3.9 Un factor fix i un aleatori amb lmer(lme4) (Opcional)
Per defecte lmer() fa servir el model no restringit. No és possible definir lmer() amb el model restringit.
library(lme4)
library(lmerTest)
res.lmer<-lmer(adhesio~cell+(1|experiment)+(1|cell:experiment),data=myData)
summary(res.lmer)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: adhesio ~ cell + (1 | experiment) + (1 | cell:experiment)
## Data: myData
##
## REML criterion at convergence: -12.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3140 -0.5718 -0.0934 0.5960 2.3662
##
## Random effects:
## Groups Name Variance Std.Dev.
## cell:experiment (Intercept) 0.00000 0.0000
## experiment (Intercept) 0.06712 0.2591
## Residual 0.03390 0.1841
## Number of obs: 60, groups: cell:experiment, 12; experiment, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.77933 0.15695 2.30615 4.966 0.0284 *
## cell2 -0.10733 0.06723 53.99995 -1.597 0.1162
## cell3 -0.20600 0.06723 53.99995 -3.064 0.0034 **
## cell4 0.04667 0.06723 53.99995 0.694 0.4906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cell2 cell3
## cell2 -0.214
## cell3 -0.214 0.500
## cell4 -0.214 0.500 0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')## 2.5 % 97.5 %
## .sig01 0.00000000 0.07439780
## .sig02 0.10643502 0.63124876
## .sigma 0.15073129 0.21789385
## (Intercept) 0.42904003 1.12962474
## cell2 -0.23781143 0.02314497
## cell3 -0.33647809 -0.07552170
## cell4 -0.08381143 0.17714497| Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| cell | 0.5753467 | 0.1917822 | 3 | 53.99995 | 5.658015 | 0.001911 |
| npar | logLik | AIC | LRT | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|
| 7 | 6.1845980 | 1.630804 | NA | NA | NA | |
| (1 | experiment) | 6 | -0.7499014 | 13.499803 | 13.869 | 1 | 0.000196 |
| (1 | cell:experiment) | 6 | 6.1845980 | -0.369196 | 0.000 | 1 | 1.000000 |
1.3.9.1 Comparacions múltiples
Per les comparacions múltiples després d’un ajust amb lmer tenim diverses alternatives, presentarem la resolució amb el package emmeans
## contrast estimate SE df t.ratio p.value
## cell1 - cell2 0.1073 0.0672 6 1.597 0.4456
## cell1 - cell3 0.2060 0.0672 6 3.064 0.0797
## cell1 - cell4 -0.0467 0.0672 6 -0.694 0.8958
## cell2 - cell3 0.0987 0.0672 6 1.468 0.5080
## cell2 - cell4 -0.1540 0.0672 6 -2.291 0.2021
## cell3 - cell4 -0.2527 0.0672 6 -3.758 0.0357
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
1.3.10 Un factor fix i un aleatori amb una observació per casella
En aquest cas no podem incloure la interacció en el model i el model lineal que resta és
\[ Y_{ij} = \mu + \alpha_i + B_j + \epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \]
on
\(Y_{ij}\) és l’observació \(k\) pel nivell \(i\) del factor \(A\) i \(j\) del factor \(B\),
\(\mu\) és la mitjana general,
\(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor A,
\(B_j\) és una variable aleatòria, no un paràmetre constant. Assumim que \(B_j \sim N(0,\sigma_B)\).
Assumim també la restricció \(\sum_{i=1}^a \alpha_i = 0\).
\(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(\epsilon_{ij} \sim N(0,\sigma)\).
Assumirem també que les v.a. \(B_j\), \((\alpha B)_{ij}\) i \(\epsilon_{ij}\) són v.a. independents.
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_E}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | \(F = \frac{MS_B}{MS_E}\) |
| Error | \(SS_E\) | \((a-1)(b-1)\) | \(MS_E = \frac{SS_E}{(a-1)(b-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
1.3.10.1 Estimació dels paràmetres del model. Components de la variància.
\(\hat{\mu} = \overline{Y}_{..}\)
\(\hat{\alpha_i} = \overline{Y}_{i.} - \overline{Y}_{..} \,\,\,\, (\alpha_i = \mu_i - \mu)\)
\(\hat{\sigma}^2 = MS_E\)
\(\hat{\sigma}^2_B = \frac{MS_B-MS_{E}}{a}\)
1.3.10.2 Resolució amb mixlm (Opcional)
El paràmetre unrestricted no és necessari donat que en aquest model no existeix versió restringida.
myData1<-myData[c(1,6,11,16,21,26,31,36,41,46,51,56),]
library(mixlm)
myData1.lm<-lm(adhesio~cell+r(experiment),unrestricted=T,data=myData1)
Anova(myData1.lm,type='III')## Analysis of variance (unrestricted model)
## Response: adhesio
## Mean Sq Sum Sq Df F value Pr(>F)
## cell 0.14 0.42 3 6.15 0.0292
## experiment 0.51 1.01 2 22.03 0.0017
## Residuals 0.02 0.14 6 - -
##
## Err.term(s) Err.df VC(SS)
## 1 cell (3) 6 fixed
## 2 experiment (3) 6 0.1206
## 3 Residuals - - 0.0229
## (VC = variance component)
##
## Expected mean squares
## cell (3) + 3 Q[1]
## experiment (3) + 4 (2)
## Residuals (3)1.4 Dissenys amb blocs
A tot experiment pot haver un factor de soroll que pugui afectar la variabilitat dels resultats però en el que no estem directament interessats.
A vegades aquest factor és desconegut i no controlable, llavors l’aleatorització és l’eina que ens permet mitigar l’efecte d’aquest factor que queda incorporat a l’error experimental. De fet les pròpies unitats experimentals són una font de variabilitat.
Altres vegades és un factor conegut i controlable llavors podem incorporar-lo al disseny i eliminar de l’error experimental la variabilitat deguda a ell. Direm que tenim un factor bloc. No estem interessats de forma específica en contrastar el factor bloc, la raó d’incorporar-lo al disseny és eliminar la variabilitat atribuïda a ell de la variabilitat residual i incrementar la potència del test sobre el factor principal.
Construir un disseny tenint en compte possibles efectes bloc rep el nom en anglés de blocking i és un dels principis formulats per Ronald Fisher juntament amb l’aleatorització i la replicació.
“Block what you can; randomize what you cannot.”
- Habitualment no són de interès les interaccions en que participa un factor bloc.
1.4.1 Disseny en blocs aleatoritzats (RCBD)
- Un dels dissenys més utilitzats en aquest cas és el disseny en blocs aleatoritzats: S’anomena disseny en blocs aleatoritzats un experiment en el qual es pren el mateix nombre de mesures per tractament dins de cada bloc i l’ordre de les mesures dins de cada bloc es decideix aleatòriament.
1.4.2 Exemple 5: Influència de l’exercici físic
Es considera que amb l’exercici físic el pH de la sang tendeix a baixar, a causa de l’acumulació de \(CO_2\) i d’ions \(H^+\) a la musculatura. Que el pH sanguini baixi per sota de 6,8 (o que pugi per sobre de 7,4) pot provocar una crisi greu, amb possibilitat de mort. Sortosament hi ha mecanismes de regulació que tendeixen a mantenir el pH en nivells adequats. En un estudi als Estats Units es volia estudiar la relació entre el pH mitjà de la sang i l’exercici físic. Es consideraren tres règims d’exercici intens A, B i C, i per qüestions administratives, els individus del estudi foren triats al atzar dintre de 6 hospitals triats al atzar. Els pH sanguinis mesurats a l’experiment foren els següents:
| Règim exercici | Hosp 1 | Hosp 2 | Hosp 3 | Hosp 4 | Hosp 5 | Hosp 6 |
|---|---|---|---|---|---|---|
| A | 7,29 | 7,20 | 7,15 | 7,23 | 7,2 | 7,1 |
| B | 7,26 | 7,16 | 7,08 | 7,11 | 7,16 | 6,89 |
| C | 7,18 | 7,13 | 6,99 | 7,03 | 7,11 | 6,90 |
pH <- c( 7.29,7.20,7.15,7.23,7.20,7.10,7.26,7.16,7.08,
7.11,7.16,6.89,7.18,7.13,6.99,7.03,7.11,6.90)
exer <- factor(c(rep(1,6),rep(2,6),rep(3,6)))
hosp <- factor(c(rep(1:6,3)))
pHdata <- data.frame(exer,hosp,pH)
boxplot(pH~exer,pHdata,col='orange')
1.4.3 RCBD: Model
1.4.3.1 Bloc Fix
\[ Y_{ij} = \mu + \alpha_i + \beta_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \]
\(\mu\) és la mitjana general,
\(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor tractament,
\(\beta_j\) és un paràmetre del model que mesura l’efecte del nivell \(j\) del factor de soroll (bloc),
\(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(E(\epsilon_{ij}) = 0\) i \(Var(\epsilon_{ij})=\sigma^2\).
La resolució és idèntica al cas de dos factors amb una rèplica per casella vist anteriorment.
\(SS_T = SS_A + SS_B + SS_E\) adaptant les fórmules al cas \(n=1\).
Tinguem en compte que l’única hipòtesi d’interès és sobre el factor primari, el factor tractament, no ens interessa el contrast sobre el bloc tot i que es podria fer.
\[ H_0: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0 \] \[ H_1: \alpha_i \ne \alpha_{i'} \mbox{ per algun } i \ne i' \,\,\, ( \mbox{ equivalentment }H_1: \exists \alpha_i \ne 0)\]
Podem utilitzar una versió simplificada de la taula ANOVA
| Font de variació | Suma de quadrats | Graus de llibertat | Quadrats mitjans | Estadístic F |
|---|---|---|---|---|
| Factor A | \(SS_A\) | \(a-1\) | \(MS_A = \frac{SS_A}{a-1}\) | \(F = \frac{MS_A}{MS_E}\) |
| Factor B | \(SS_B\) | \(b-1\) | \(MS_B = \frac{SS_B}{b-1}\) | |
| Error | \(SS_E\) | \((a-1)(b-1)\) | \(MS_E = \frac{SS_E}{(a-1)(b-1)}\) | |
| Total | \(SS_T\) | \(N-1\) |
1.4.3.2 Bloc aleatori
És molt habitual que tal i com passa a l’exemple 5 el factor bloc sigui un factor aleatori. En aquest cas la taula i les consideracions són idèntiques. Tal sols cal tenir en compte que si es vol contrastar l’efecte del bloc cal interpretar de forma adequada la hipòtesis nul·la corresponent i adaptar els paràmetres del model.
\[ Y_{ij} = \mu + \alpha_i + B_j +\epsilon_{ij} \,\,\,\,\,\, i=1,...,a; \,\,\,\, j=1,...,b \]
\(\mu\) és la mitjana general,
\(\alpha_i\) és un paràmetre del model que mesura l’efecte del nivell \(i\) del factor tractament,
\(B_j\) és una variable aleatòria \(B_j \sim N(0,\sigma_B)\) que mesura la variabilitat aportada pel factor de soroll (bloc),
\(\epsilon_{ij}\) és l’error aleatori corresponent a l’observació \(i,j\) i dels quals suposem que són independents i idènticament distribuïts amb \(E(\epsilon_{ij}) = 0\) i \(Var(\epsilon_{ij})=\sigma^2\). Són independents de les \(B_j\).
La resolució és idèntica perquè no ens interessa contrastar el factor bloc.
1.4.4 Resolució de l’exemple 5
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| exer | 2 | 0.0584111 | 0.0292056 | 18.49754 | 0.0004363 |
| hosp | 5 | 0.1354944 | 0.0270989 | 17.16327 | 0.0001267 |
| Residuals | 10 | 0.0157889 | 0.0015789 | NA | NA |
##
## Study: pHdata.aov ~ "exer"
##
## HSD Test for pH
##
## Mean Square Error: 0.001578889
##
## exer, means
##
## pH std r se Min Max Q25 Q50 Q75
## 1 7.195000 0.06534524 6 0.01622184 7.10 7.29 7.1625 7.200 7.2225
## 2 7.110000 0.12393547 6 0.01622184 6.89 7.26 7.0875 7.135 7.1600
## 3 7.056667 0.10308572 6 0.01622184 6.90 7.18 7.0000 7.070 7.1250
##
## Alpha: 0.05 ; DF Error: 10
## Critical Value of Studentized Range: 3.876777
##
## Minimun Significant Difference: 0.06288846
##
## Treatments with the same letter are not significantly different.
##
## pH groups
## 1 7.195000 a
## 2 7.110000 b
## 3 7.056667 bAnem a veure que succeeix si no incloem el factor hospital
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| exer | 2 | 0.0584111 | 0.0292056 | 2.89578 | 0.0864035 |
| Residuals | 15 | 0.1512833 | 0.0100856 | NA | NA |
1.4.4.1 Consideracions
El model assumeix que no hi ha interacció. Si la suposició és falsa és perillós de cara a interpretar els resultats. veure la consideració anterior sobre la prova d’aditivitat de Tukey.
Es un disseny molt eficient respecte al número de casos necessaris.
Són vàlids tots els mètodes de comparacions múltiples.
L’estimació de paràmetres és idèntica a la vista anteriorment.
Cal verificar les suposicions per poder aplicar l’ANOVA.
1.4.5 Resolució de l’exemple 5 amb els packages nlme i lme4
1.4.5.1 Package nlme
## Linear mixed-effects model fit by REML
## Data: pHdata
## AIC BIC logLik
## -24.60822 -21.06797 17.30411
##
## Random effects:
## Formula: ~1 | hosp
## (Intercept) Residual
## StdDev: 0.09223159 0.03973524
##
## Fixed effects: pH ~ exer
## Value Std.Error DF t-value p-value
## (Intercept) 7.195000 0.04099910 10 175.49167 0.0000
## exer2 -0.085000 0.02294115 10 -3.70513 0.0041
## exer3 -0.138333 0.02294115 10 -6.02992 0.0001
## Correlation:
## (Intr) exer2
## exer2 -0.28
## exer3 -0.28 0.50
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -1.81043774 -0.54870477 0.09034733 0.46540220 1.33538484
##
## Number of Observations: 18
## Number of Groups: 6## Approximate 95% confidence intervals
##
## Fixed effects:
## lower est. upper
## (Intercept) 7.1036483 7.1950000 7.28635168
## exer2 -0.1361161 -0.0850000 -0.03388394
## exer3 -0.1894494 -0.1383333 -0.08721727
##
## Random Effects:
## Level: hosp
## lower est. upper
## sd((Intercept)) 0.04773274 0.09223159 0.1782145
##
## Within-group standard error:
## lower est. upper
## 0.02563547 0.03973524 0.06159001| numDF | denDF | F-value | p-value | |
|---|---|---|---|---|
| (Intercept) | 1 | 10 | 33678.19389 | 0.0000000 |
| exer | 2 | 10 | 18.49754 | 0.0004363 |
1.4.5.2 Package lme4
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: pH ~ exer + (1 | hosp)
## Data: pHdata
##
## REML criterion at convergence: -34.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.81044 -0.54870 0.09035 0.46540 1.33538
##
## Random effects:
## Groups Name Variance Std.Dev.
## hosp (Intercept) 0.008507 0.09223
## Residual 0.001579 0.03974
## Number of obs: 18, groups: hosp, 6
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 7.19500 0.04100 6.19114 175.492 1.11e-12 ***
## exer2 -0.08500 0.02294 10.00000 -3.705 0.004074 **
## exer3 -0.13833 0.02294 10.00000 -6.030 0.000127 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) exer2
## exer2 -0.280
## exer3 -0.280 0.500## 2.5 % 97.5 %
## .sig01 0.04932455 0.17271691
## .sigma 0.02548100 0.05753568
## (Intercept) 7.11029296 7.27970703
## exer2 -0.12956104 -0.04043896
## exer3 -0.18289437 -0.09377229| Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| exer | 0.0584111 | 0.0292056 | 2 | 10 | 18.49754 | 0.0004363 |
| npar | logLik | AIC | LRT | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|
| 5 | 17.304109 | -24.60822 | NA | NA | NA | |
| (1 | hosp) | 4 | 9.404553 | -10.80911 | 15.79911 | 1 | 7.04e-05 |
1.4.6 Resolució de l’exemple 4 com un model on experiment és un bloc
Considerearem que el factor experiment de l’exemple 4 és un bloc i no inclourem la interacció amb el factor fix cell.
1.4.6.1 Resolució amb aov
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| cell | 3 | 0.5753467 | 0.1917822 | 5.658024 | 0.001911 |
| experiment | 2 | 2.7526633 | 1.3763317 | 40.605004 | 0.000000 |
| Residuals | 54 | 1.8303633 | 0.0338956 | NA | NA |
##
## Study: ex4.bloc ~ "cell"
##
## HSD Test for adhesio
##
## Mean Square Error: 0.03389562
##
## cell, means
##
## adhesio std r se Min Max Q25 Q50 Q75
## 1 0.7793333 0.2731684 15 0.04753638 0.25 1.31 0.690 0.77 0.940
## 2 0.6720000 0.3114299 15 0.04753638 0.18 1.24 0.420 0.65 0.900
## 3 0.5733333 0.2619887 15 0.04753638 0.15 0.94 0.360 0.69 0.775
## 4 0.8260000 0.2951465 15 0.04753638 0.19 1.36 0.615 0.84 1.025
##
## Alpha: 0.05 ; DF Error: 54
## Critical Value of Studentized Range: 3.748904
##
## Comparison between treatments means
##
## difference pvalue signif. LCL UCL
## 1 - 2 0.10733333 0.3891 -0.07087601 0.28554267
## 1 - 3 0.20600000 0.0174 * 0.02779066 0.38420934
## 1 - 4 -0.04666667 0.8989 -0.22487601 0.13154267
## 2 - 3 0.09866667 0.4638 -0.07954267 0.27687601
## 2 - 4 -0.15400000 0.1129 -0.33220934 0.02420934
## 3 - 4 -0.25266667 0.0023 ** -0.43087601 -0.074457331.4.6.2 Resolució amb el package mixlm
## Analysis of variance (unrestricted model)
## Response: adhesio
## Mean Sq Sum Sq Df F value Pr(>F)
## cell 0.19 0.58 3 5.66 0.0019
## experiment 1.38 2.75 2 40.61 0.0000
## Residuals 0.03 1.83 54 - -
##
## Err.term(s) Err.df VC(SS)
## 1 cell (3) 54 fixed
## 2 experiment (3) 54 0.0671
## 3 Residuals - - 0.0339
## (VC = variance component)
##
## Expected mean squares
## cell (3) + 15 Q[1]
## experiment (3) + 20 (2)
## Residuals (3)1.4.6.3 Resolució amb el package lme4
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: adhesio ~ cell + (1 | experiment)
## Data: myData
##
## REML criterion at convergence: -12.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3140 -0.5718 -0.0934 0.5960 2.3662
##
## Random effects:
## Groups Name Variance Std.Dev.
## experiment (Intercept) 0.06712 0.2591
## Residual 0.03390 0.1841
## Number of obs: 60, groups: experiment, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.77933 0.15695 2.30598 4.965 0.0284 *
## cell2 -0.10733 0.06723 54.00000 -1.597 0.1162
## cell3 -0.20600 0.06723 54.00000 -3.064 0.0034 **
## cell4 0.04667 0.06723 54.00000 0.694 0.4906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) cell2 cell3
## cell2 -0.214
## cell3 -0.214 0.500
## cell4 -0.214 0.500 0.500| Sum Sq | Mean Sq | NumDF | DenDF | F value | Pr(>F) | |
|---|---|---|---|---|---|---|
| cell | 0.5753467 | 0.1917822 | 3 | 54 | 5.658024 | 0.001911 |
| npar | logLik | AIC | LRT | Df | Pr(>Chisq) | |
|---|---|---|---|---|---|---|
| 6 | 6.184598 | -0.369196 | NA | NA | NA | |
| (1 | experiment) | 5 | -16.179175 | 42.358350 | 44.72755 | 1 | 0 |
## contrast estimate SE df t.ratio p.value
## cell1 - cell2 0.1073 0.0672 54 1.597 0.3891
## cell1 - cell3 0.2060 0.0672 54 3.064 0.0174
## cell1 - cell4 -0.0467 0.0672 54 -0.694 0.8989
## cell2 - cell3 0.0987 0.0672 54 1.468 0.4638
## cell2 - cell4 -0.1540 0.0672 54 -2.291 0.1129
## cell3 - cell4 -0.2527 0.0672 54 -3.758 0.0023
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 4 estimates
1.5 El cas no balancejat
Eliminarem una observació del conjunt de dades del primer exemple i veurem què passa al fer un Anova intercanviant l’ordre dels factors.
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| vacuna | 3 | 141.498359 | 47.1661197 | 11.3616797 | 0.0000236 |
| dosis | 2 | 8.178113 | 4.0890564 | 0.9849983 | 0.3835536 |
| vacuna:dosis | 6 | 5.375372 | 0.8958953 | 0.2158091 | 0.9692567 |
| Residuals | 35 | 145.296667 | 4.1513333 | NA | NA |
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| dosis | 2 | 7.902302 | 3.9511512 | 0.9517788 | 0.3958190 |
| vacuna | 3 | 141.774170 | 47.2580565 | 11.3838261 | 0.0000232 |
| dosis:vacuna | 6 | 5.375372 | 0.8958953 | 0.2158091 | 0.9692567 |
| Residuals | 35 | 145.296667 | 4.1513333 | NA | NA |
Per què obtenim diferències en les taules Anova?
1.5.1 Dissenys balancejats vs no balancejats
Hi han moltes raons per preferir un disseny balancejat:
Els contrasts són molt més robusts a la suposició de normalitat i homoscedasticitat quan les mides mostrals són iguals.
L’estimació dels components de la variància és molt més complicada amb mides mostrals diferents.
Però més enllà d’aquestes raons en un disseny no balancejat no hi ha una única partició additiva de la suma de quadrats total. Això implica que hi han diferents formes de calcular les sumes de quadrats. En general s’anomenen Tipus I, II i III.
R per defecte utilitza el tipus I, a diferència d’altres programes com per exemple SPSS que utilitza el tipus III.
Fem servir la notació SS(A,B,AB) per indicar un model amb els factors A, B i la seva interacció AB; SS(A,B) el model sense la interacció; SS(B,AB) el model sense els efectes del factor A i així successivament. Llavors els diferents tipus de sumes de quadrats assignen els següents valors com a sumes de quadrats:
Tipus I (seqüencial) model A*B
- SS(A) pel factor A
- SS(B|A) = SS(A,B)-SS(A) pel factor B
- SS(AB|B,A) = SS(A,B,AB)-SS(A,B) per a la interacció AB
Els contrasts que es fan són: sobre el factor A, sobre el factor B després d’haver controlat pel factor A i sobre la interacció després d’haver controlat pels dos efectes principals.
Tipus II
- SS(A|B) = SS(A,B)-SS(B) per el factor A
- SS(B|A) = SS(A,B)-SS(A) per el factor B
Tipus III
- SS(A|B,AB) = SS(A,B,AB)-SS(B,AB) per el factor A
- SS(B|A,AB) = SS(A,B,AB)-SS(A,AB) per el factor B
- SS(AB|B,A) = SS(A,B,AB)-SS(A,B) per a la interacció AB
La suma de quadrats de tipus I utilitza un procediment jeràrquic pel càlcul de les sumes de quadrats i l’ordre dels termes en el model té importància.
Molts autors recomanen utilitzar el tipus III que no està influenciat per l’ordre dels termes.
1.5.2 Comprovació de les fòrmules de les sumes de quadrats
Calculem manualment les sumes de quadrats dels diferents models
mod0<-lm(anticos~1,immuno_1)
SST<-sum(mod0$residuals^2)
modA<-lm(anticos~vacuna,immuno_1)
modB<-lm(anticos~dosis,immuno_1)
SSRA<-sum(modA$residuals^2)
SSA<-SST-SSRA
SSRB<-sum(modB$residuals^2)
SSB<-SST-SSRB
modA_B<-lm(anticos~vacuna+dosis,immuno_1)
SSRA_B<-sum(modA_B$residuals^2)
SSA_B<-SST-SSRA_B
modAB<-lm(anticos~vacuna*dosis,immuno_1)
SSRAB<-sum(modAB$residuals^2)
SSAB<-SST-SSRAB
X<-model.matrix(lm(anticos~vacuna*dosis,immuno_1))
interi<-X[,7:12]
modnB<-lm(anticos~vacuna+interi,immuno_1)
modnA<-lm(anticos~dosis+interi,immuno_1)
SSRnA<-sum(modnA$residuals^2)
SSRnB<-sum(modnB$residuals^2)
SSnA<-SST-SSRnA
SSnB<-SST-SSRnB
c(SST,SSA,SSB,SSA_B,SSAB,SSnA,SSnB)## [1] 300.348511 141.498359 7.902302 149.676472 155.051844 12.099152 146.768099Segons les taules anteriors les sumes de quadrats han de ser:
1.5.2.1 Tipus I amb vacuna*dosi
- SS(vacuna): SSA = 141.4983591
- SS(dosi): SSA_B - SSA = 8.1781128
- SS(interacció): SSAB-SSA_B = 5.375372
1.5.2.2 Tipus I amb dosi*vacuna
- SS(dosi): SSB = 7.9023023
- SS(vacuna): SSA_B - SSB = 141.7741696
- SS(interacció): SSAB - SSA_B =5.375372
1.5.2.3 Tipus II sense interacció
- SS(vacuna): SSA_B - SSB = 141.7741696
- SS(dosi): SSA_B - SSA = 8.1781128
1.5.2.4 Tipus III
- SS(vacuna): SSAB - SSnA = 142.9526923
- SS(dosi): SSAB - SSnB = 8.2837445
- SS(interacció): SSAB - SSA-B = 5.375372
1.5.3 Suma de quadrats tipus III amb R
1.5.3.1 Amb el package sasLM
## Response : anticos
## Df Sum Sq Mean Sq F value Pr(>F)
## MODEL 11 155.052 14.096 3.3954 0.002813 **
## dosis 2 8.284 4.142 0.9977 0.378963
## vacuna 3 142.953 47.651 11.4785 2.167e-05 ***
## dosis:vacuna 6 5.375 0.896 0.2158 0.969257
## RESIDUALS 35 145.297 4.151
## CORRECTED TOTAL 46 300.349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Response : anticos
## Df Sum Sq Mean Sq F value Pr(>F)
## MODEL 11 155.052 14.096 3.3954 0.002813 **
## vacuna 3 142.953 47.651 11.4785 2.167e-05 ***
## dosis 2 8.284 4.142 0.9977 0.378963
## vacuna:dosis 6 5.375 0.896 0.2158 0.969257
## RESIDUALS 35 145.297 4.151
## CORRECTED TOTAL 46 300.349
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 11.5.3.2 Amb el package car
## Funció Anova del package car
immuno_1.aov<-aov(anticos~vacuna*dosis, data=immuno_1,contrasts=list(vacuna=contr.sum,dosis=contr.sum))
car::Anova(immuno_1.aov,type=3)| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 1170.163604 | 1 | 281.8765706 | 0.0000000 |
| vacuna | 142.952692 | 3 | 11.4784561 | 0.0000217 |
| dosis | 8.283744 | 2 | 0.9977210 | 0.3789631 |
| vacuna:dosis | 5.375372 | 6 | 0.2158091 | 0.9692567 |
| Residuals | 145.296667 | 35 | NA | NA |
immuno_1.aov<-aov(anticos~dosis*vacuna, data=immuno_1,contrasts=list(vacuna=contr.sum,dosis=contr.sum))
car::Anova(immuno_1.aov,type=3)| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 1170.163604 | 1 | 281.8765706 | 0.0000000 |
| dosis | 8.283744 | 2 | 0.9977210 | 0.3789631 |
| vacuna | 142.952692 | 3 | 11.4784561 | 0.0000217 |
| dosis:vacuna | 5.375372 | 6 | 0.2158091 | 0.9692567 |
| Residuals | 145.296667 | 35 | NA | NA |
Pots trobar més informació en el següent enllaç
1.6 Dissenys amb més de dos factors
La teoria dels models lineals pot ser estesa a models amb més de dos factors, per exemple un model amb 3 factors fixos seria el següent
\[ Y_{ijkl} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha \beta )_{ij}+ (\alpha \gamma )_{ik}+(\beta \gamma )_{jk}+(\alpha \beta \gamma)_{ijk}+\epsilon_{ijkl} \\ i=1,...,a; \,\,\,\, j=1,...,b; \,\,\,\, k=1,...,c;\,\,\,\, l=1,...,n_{ijk} \]
amb
\[ \sum_{i=1}^a \alpha_i = 0 \,\, , \,\,\sum_{j=1}^b \beta_j = 0 \,\, , \,\,\sum_{k=1}^c \gamma_k = 0 \]
\[ \sum_{i=1}^a (\alpha \beta )_{ij} = 0 \,\, , \,\, \sum_{j=1}^b (\alpha \beta )_{ij} = 0 \,\, , \,\, \sum_{k=1}^c (\alpha \gamma )_{ik} = 0 \,\, , \mbox{ etc. } \]
\[ \sum_{i=1}^a \sum_{j=1}^b (\alpha \beta \gamma )_{ijk} = 0 \,\, , \,\, \mbox{ etc. } \]
Si dos o més dels factors son aleatoris per a alguns termes no hi haurà un quocient F adequat per contrastar la seva significació. El problema es torna quasi intractable si tenim múltiples factors aleatoris i un disseny no balancejat.
Com a mostra veiem la següent taula extreta del llibre de G. Quinn i M. Keough Experimental Design and Data Analysis for Biologists, Cambridge University Press, 2002.
A la taula \(\sigma^2_\alpha\) es refereix a la variància aportada si el factor A és aleatori o bé a \(\frac{\sum_{i=1}^a \alpha^2}{a-1}\) si el factor és fix.
1.6.1 Exemple i resolució amb R
colesterol<-read.table('../dades/cholesterol.csv', sep=';',dec=',',header=T,stringsAsFactors = T)
colesterol.aov<-aov(cholesterol~drug*sex*risk,data=colesterol)
anova(colesterol.aov)| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| drug | 2 | 1.3952111 | 0.6976056 | 5.4661048 | 0.0065953 |
| sex | 1 | 0.0032000 | 0.0032000 | 0.0250737 | 0.8747154 |
| risk | 1 | 0.0000889 | 0.0000889 | 0.0006965 | 0.9790330 |
| drug:sex | 2 | 0.0079000 | 0.0039500 | 0.0309503 | 0.9695392 |
| drug:risk | 2 | 0.7256778 | 0.3628389 | 2.8430327 | 0.0661211 |
| sex:risk | 1 | 0.0060500 | 0.0060500 | 0.0474049 | 0.8283808 |
| drug:sex:risk | 2 | 0.0399000 | 0.0199500 | 0.1563187 | 0.8556338 |
| Residuals | 60 | 7.6574333 | 0.1276239 | NA | NA |